DATA_SECTION
  init_int styr  //first year of model (i.e., model goes from styr.0000000000000000001)
  init_int endyr //last year of model (i.e., model goes to endyr.999999999999999999)
  int n_model_years //total number of years in model
  !! n_model_years = endyr-styr+1;
  init_int n_intrvls_per_yr //how many intervals in each year
  init_matrix fishery_sample_prob(1,n_model_years,1,n_intrvls_per_yr)  //sample probabilities for fishery in each interval
  init_matrix survey_sample_prob(1,n_model_years,1,n_intrvls_per_yr)   //sample probabilities for survey in each interval
  init_matrix survey_effort(1,n_model_years,1,n_intrvls_per_yr)        //effort for survey in each interval
  vector fishery_sample_prob_vec(1,n_model_years*n_intrvls_per_yr)     //sample probabilities for fishery in each interval, vector form
  vector survey_sample_prob_vec(1,n_model_years*n_intrvls_per_yr)      //sample probabilities for survey in each interval, vector form
  vector survey_effort_vec(1,n_model_years*n_intrvls_per_yr)           //effort for survey in each interval, vector form
  init_int max_age  //maximum age for each cohort considered in model
  int first_cohort //first cohort has observation for first model year, so it is age=(max_age-1).0 at initial model time
  int last_cohort //last cohort has observation for last model year, so it is age=1.0 at end of model time
  !! first_cohort = styr-max_age+1;  
  !! last_cohort = endyr; 


  init_int fishery_n_obs  //number of fishery observations
  init_matrix fishery_dat(1,fishery_n_obs,1,4) //all fishery data

  init_int survey_n_obs  //number of survey observations
  init_matrix survey_dat(1,survey_n_obs,1,4) //all survey data

  //fishery data parts
  vector fishdat_tend(1,fishery_n_obs)			//ends of intervals
  vector fishdat_age_mid(1,fishery_n_obs)		//ages at middle of intervals
  vector fishdat_catch(1,fishery_n_obs)			//numbers caught in intervals
  vector fishdat_cohort(1,fishery_n_obs)		//cohort id for each obs
  ivector fishdat_cohort_index(1,fishery_n_obs)	//cohort index for each obs  (row in ragged arrays) 
  ivector fishdat_cohort_intrvl_index(1,fishery_n_obs)	//cohort intrvl index for each obs (column in ragged arrays)
  ivector fishdat_year(1,fishery_n_obs)			//calendar year id for each obs
  ivector fishdat_month(1,fishery_n_obs)		//calendar month (1-12) id  for each obs

  //survey data parts
  vector survey_tend(1,survey_n_obs)			//ends of intervals                     
  vector survey_age_mid(1,survey_n_obs)			//ages at middle of intervals           
  vector survey_catch(1,survey_n_obs)			//numbers caught in intervals           
  vector survey_cohort(1,survey_n_obs)			//cohort id for each obs                
  ivector survey_cohort_index(1,survey_n_obs)	//cohort index for each obs
  ivector survey_cohort_intrvl_index(1,fishery_n_obs)	//cohort intrvl index for each obs (column in ragged arrays)
  ivector survey_year(1,survey_n_obs)			//calendar year id for each obs         
  ivector survey_month(1,survey_n_obs)			//calendar month (1-12) id  for each obs

  int n_cohorts //total number of cohorts
  !! n_cohorts = n_model_years+max_age-1; 
  !! cout << n_cohorts << " " << first_cohort << " " << last_cohort << " " << last_cohort - first_cohort + 1 << endl;
  number Min //natural mortality
  ivector   n_intrvls_cohort(1,n_cohorts)  //total number of model intervals for each cohort
  ivector beg_intrvl_cohort(1,n_cohorts)  //first model interval for each cohort
  ivector end_intrvl_cohort(1,n_cohorts)  //last model interval for each cohort
  int n_model_intrvls // total number of intervals in the model (nyrs X 12 if monthly)
  //int minK
  //int maxK
  //int nages
  //number offset

  LOCAL_CALCS
  n_model_intrvls = (endyr-styr+1)*n_intrvls_per_yr;
  // Get data into vectors (from matrix
  fishdat_tend    = column(fishery_dat,1);
  fishdat_age_mid = column(fishery_dat,2);
  fishdat_catch     = column(fishery_dat,3);
  fishdat_cohort     = column(fishery_dat,4);
  for (int i=1;i<=fishery_n_obs;i++)
  {
    fishdat_year(i) = int(fishdat_tend(i) - .001);    
    fishdat_month(i) = n_intrvls_per_yr*(fishdat_tend(i) - fishdat_year(i)+.001);    
  }
  survey_sample_prob_vec.initialize();
  fishery_sample_prob_vec.initialize();
  survey_effort_vec.initialize();
  for (int i=1;i<=n_model_years;i++)
  {
  /*
   	cout << n_model_years << endl;
  	cout << i <<endl;
  	cout << (i-1)*n_intrvls_per_yr+1 << " " << i*n_intrvls_per_yr << endl;
  	cout << survey_sample_prob(i) << endl;
  	cout << survey_sample_prob_vec((i-1)*n_intrvls_per_yr+1,i*n_intrvls_per_yr) << endl;
  */
  	for(int j = 1; j <=n_intrvls_per_yr;j++)
  	{
      survey_sample_prob_vec((i-1)*n_intrvls_per_yr+j) = survey_sample_prob(i,j);
  	  //cout << survey_sample_prob_vec((i-1)*n_intrvls_per_yr+j) << endl;
      survey_effort_vec((i-1)*n_intrvls_per_yr+j) = survey_effort(i,j);
      fishery_sample_prob_vec((i-1)*n_intrvls_per_yr+j) = fishery_sample_prob(i,j);
    }
  }

  survey_tend    = column(survey_dat,1);
  survey_age_mid = column(survey_dat,2);
  survey_catch     = column(survey_dat,3);
  survey_cohort     = column(survey_dat,4);
  for (int i=1;i<=survey_n_obs;i++)
  {
    survey_year(i) = int(survey_tend(i) - .001);    
    survey_month(i) = n_intrvls_per_yr*(survey_tend(i) - survey_year(i)+.001);    
  }
  // cout <<fishdat_year<<endl;
  //cout <<survey_year<<endl;

  //cout <<fishdat_month<<endl;
  //cout <<survey_month<<endl;
  // Vector of interval definitions
  
  // Compute the number of observations by 
  //nages = max_age;
  //maxK  = 1+endyr-styr+max_age;
  //minK  = 1;
  Min=0.2;
  //offset = 0.0;
  end_intrvl_cohort.initialize();
  beg_intrvl_cohort.initialize();
  n_intrvls_cohort.initialize();
  for(int k=1;k<=max_age;k++)
  {
    beg_intrvl_cohort(k) = 1;
    end_intrvl_cohort(k) = k*n_intrvls_per_yr;
    n_intrvls_cohort(k)  = end_intrvl_cohort(k) - beg_intrvl_cohort(k) + 1;
  }
  for(int k=max_age;k<=n_cohorts;k++)
  {
    beg_intrvl_cohort(k) = (k-max_age)*n_intrvls_per_yr + 1 ;
    if(max_age*n_intrvls_per_yr < n_model_intrvls - beg_intrvl_cohort(k)+1.)
      end_intrvl_cohort(k) = beg_intrvl_cohort(k) + max_age*n_intrvls_per_yr;
    else
      end_intrvl_cohort(k) = n_model_intrvls;
    // Number of intervals for each cohort
    n_intrvls_cohort(k)  = end_intrvl_cohort(k) - beg_intrvl_cohort(k) + 1;
  } 
  // cout<<n_intrvls_cohort<<endl;
  // cout<<beg_intrvl_cohort<<endl;
  // cout<<end_intrvl_cohort<<endl;
  // cout<<n_cohorts.indexmax()<<endl;
 END_CALCS
  vector Ages(1,max_age*n_intrvls_per_yr);       				//all possible observed ages for every interval
  // matrix Dindex_F(1,n_cohorts,1,n_intrvls_cohort);
  // matrix Dindex_S(1,n_cohorts,1,n_intrvls_cohort);
  // imatrix fishdat_year(1,n_cohorts,1,n_intrvls_cohort)
  // imatrix survey_year(1,n_cohorts,1,n_intrvls_cohort)
  matrix ages_cohort(1,n_cohorts,1,n_intrvls_cohort); 			//for each cohort, ages in each possible cohort interval
  imatrix cal_yr_cohort(1,n_cohorts,1,n_intrvls_cohort);		//for each cohort, the calendar year associated with each possible cohort interval
  imatrix cal_intrvl_cohort(1,n_cohorts,1,n_intrvls_cohort);	//for each cohort, the calendar interval associated with each possible cohort interval
  matrix fishery_catch_cohort(1,n_cohorts,1,n_intrvls_cohort); 	//for each cohort, catch in fishery in each possible cohort interval
  matrix survey_catch_cohort(1,n_cohorts,1,n_intrvls_cohort);	//for each cohort, catch in survey in each possible cohort interval
  vector total_fishery_catch_cohort(1,n_cohorts);				//for each cohort, the total amount caught in the fishery
  vector total_survey_catch_cohort(1,n_cohorts);				//for each cohort, the total amount caught in the survey

 LOCAL_CALCS
  for (int i=1;i<=n_intrvls_per_yr*max_age;i++)
    Ages(i) = double((i-.5)/double(n_intrvls_per_yr));
  // cout << Ages<<endl;exit(1);
  fishdat_year.initialize();
  survey_year.initialize();
  total_survey_catch_cohort.initialize();
  survey_catch_cohort.initialize();
  total_fishery_catch_cohort.initialize();
  fishery_catch_cohort.initialize();
  int ii;
  for (int i=1;i<=survey_n_obs;i++)
  {
    //int k  = survey_cohort(i)-styr+max_age;
   //cout <<"i: " << i << " of " << survey_n_obs <<endl;
   survey_cohort_index(i) = survey_cohort(i)-styr+max_age;
   if (survey_cohort_index(i) <= max_age)
      ii = n_intrvls_per_yr*(survey_tend(i) - styr +.001);
    else
      ii = n_intrvls_per_yr*(survey_tend(i) - survey_cohort(i)+.001);

    // Dindex_S(k,ii) = i;
    survey_cohort_intrvl_index(i) = ii;
    survey_catch_cohort(survey_cohort_index(i),ii) += survey_catch(i);
    total_survey_catch_cohort(survey_cohort_index(i))    += survey_catch(i);
  }
  for (int i=1;i<=fishery_n_obs;i++)
  {
    //int k  = fishdat_cohort(i)-styr+max_age;
    fishdat_cohort_index(i) = fishdat_cohort(i)-styr+max_age;
    if (fishdat_cohort_index(i) <= max_age)
      ii = n_intrvls_per_yr*(fishdat_tend(i) - styr +.001);
    else
      ii = n_intrvls_per_yr*(fishdat_tend(i) - fishdat_cohort(i)+.001);
    fishdat_cohort_intrvl_index(i) = ii;
    fishery_catch_cohort(fishdat_cohort_index(i),ii) += fishdat_catch(i);
    total_fishery_catch_cohort(fishdat_cohort_index(i)) += fishdat_catch(i);
  }
  for (int k=1;k<=n_cohorts;k++)
  {
   for (int i=1;i<=n_intrvls_cohort(k);i++)
    {
      if (k<=max_age)
      {
        // Indexes age for the cohort by interval
        ages_cohort(k,i) = n_intrvls_per_yr*(max_age - k) + i; 
        cal_yr_cohort(k,i) = styr + int((double(i)-.5)/double(n_intrvls_per_yr));
      }
      else
      {
        ages_cohort(k,i) =  i ;
        cal_yr_cohort(k,i) = k-max_age +styr + int((double(i)-.5)/double(n_intrvls_per_yr));
      }
      cal_intrvl_cohort(k,i) = n_intrvls_per_yr*(cal_yr_cohort(k,1) - styr) + i;
    }
  }
  // Print out the Catch by cohort
  // cout << cal_intrvl_cohort <<endl;exit(1);
  //cout <<"Hear"<<endl;
  //exit(1);
 END_CALCS
  init_int N0_phase
  init_int F_phase
  init_int q_phase
  init_int fishery_sel_phase
  init_int survey_sel_phase
  init_int M_phase

INITIALIZATION_SECTION
  log_N0 13.52393
  log_F -1.203973
  fishery_sel50 5.
  fishery_selslope 1.
  survey_sel50 2.
  survey_selslope 5.
  log_q -4.199705
  log_M -1.609438

PARAMETER_SECTION
  init_vector log_N0(1,n_cohorts,N0_phase) //initial abundance for each cohort
  init_vector log_F(1,n_model_years,F_phase)      //fishing mortality for each calendar year
  init_number log_q(q_phase)              //survey catchability

  matrix F_fishery_cohort(1,n_cohorts,1,n_intrvls_cohort) //fishing mortality for each interval for each cohort (including selectivity effects)
  matrix F_survey_cohort(1,n_cohorts,1,n_intrvls_cohort) //survey mortality for each interval for each cohort (including selectivity effects)
  matrix Z_cohort(1,n_cohorts,1,n_intrvls_cohort)	  //Z (F_f + F_s + M) for each interval for each cohort
  matrix survival_cohort(1,n_cohorts,1,n_intrvls_cohort) //survival for each interval for each cohort (exp(-Z/n_intervals_per_yr))
  matrix cumulative_survival_cohort(1,n_cohorts,1,n_intrvls_cohort) //cumulative survival for each interval for each cohort
  vector cumulative_nat_mort_cohort(1,n_cohorts) //cumulative natural mortality for each cohort
  vector F(1,n_model_years)
  number q
  init_number fishery_sel50(fishery_sel_phase)		//fishery logistic selectivity 50% selected parameter
  init_number fishery_selslope(fishery_sel_phase)	//fishery logistic selectivity slope parameter
  init_number survey_sel50(survey_sel_phase)		//survey logistic selectivity 50% selected parameter
  init_number survey_selslope(survey_sel_phase)		//survey logistic selectivity slope parameter
  vector fishery_select(1,max_age*n_intrvls_per_yr) //fishery selectivity for each possible observed age
  vector survey_select(1,max_age*n_intrvls_per_yr)	//survey selectivity for each possible observed age 
  // matrix S_trans(1,ntyp+1,1,ntyp+1) // Survival transition matrix from M and each fishery (or survey)
  vector N0(1,n_cohorts)	//Abundance at time 0 for each cohort
  init_number log_M(M_phase)
  number M
  init_number log_init_N0(2) //the average recruitment that yields the observations from cohorts greater than age zero 
 							 //at the beginning of the model
  init_number log_F0(-1)     //the average fully selected fishing mortality that yields the observations from cohorts 
  							 //greater than age zero at the beginning of the model
  
  number log_exploit
  number ll_i
  number ll
  objective_function_value nll

PROCEDURE_SECTION
  F = mfexp(log_F);
  q = mfexp(log_q);
  M = mfexp(log_M);
  F_fishery_cohort.initialize();
  // cumulative_survival_cohort = 1.;
  // Compute selectivity (over all monthly ages?)
  for (int i=1; i<=max_age*n_intrvls_per_yr; i++)
  {
    fishery_select(i) = 1./(1+mfexp(-fishery_selslope*(Ages(i) - fishery_sel50)));
    survey_select(i) = 1./(1+mfexp(-survey_selslope*(Ages(i) - survey_sel50)));
  }
  //for(i=1;i<=N_int(max_age);i++){
  //  N0(max_age-i) = mfexp(log_init_N0)*mexp(-(M+mfexp(log_F0)*fishery_selec(
  //}
  for (int k=1; k<=n_cohorts; k++)
  {
    //int iy   = cal_yr_cohort(k,1); // assign temporary if using multiple times for overhead reasons
    int iy = cal_yr_cohort(k,1)- styr + 1;
    int i_intrvl = cal_intrvl_cohort(k,1);
    int iage = ages_cohort(k,1);
    F_fishery_cohort(k,1) = fishery_select(iage) * F(iy) * fishery_sample_prob_vec(i_intrvl);
    F_survey_cohort(k,1) = survey_select(iage) * q * survey_sample_prob_vec(i_intrvl) * survey_effort_vec(i_intrvl);
    Z_cohort(k,1) = F_fishery_cohort(k,1) + F_survey_cohort(k,1) + M;
    survival_cohort(k,1) = mfexp(-Z_cohort(k,1)/double(n_intrvls_per_yr));
    cumulative_survival_cohort(k,1) = survival_cohort(k,1);
    cumulative_nat_mort_cohort(k) = (1. - survival_cohort(k,1))*M/Z_cohort(k,1);
    for (int i=2; i<=n_intrvls_cohort(k); i++)
    {
      //iy   = cal_yr_cohort(k,i); // assign temporary if using multiple times for overhead reasons
      int iy = cal_yr_cohort(k,i)- styr + 1;
      i_intrvl = cal_intrvl_cohort(k,i);
      iage = ages_cohort(k,i);
      F_fishery_cohort(k,i) = fishery_select(iage) * F(iy) * fishery_sample_prob_vec(i_intrvl);
      F_survey_cohort(k,i) = survey_select(iage) * q * survey_sample_prob_vec(i_intrvl) * survey_effort_vec(i_intrvl);
      Z_cohort(k,i) = F_fishery_cohort(k,i) + F_survey_cohort(k,i) + M;
      survival_cohort(k,i) = mfexp(-Z_cohort(k,i)/double(n_intrvls_per_yr));
      cumulative_survival_cohort(k,i) = cumulative_survival_cohort(k,i-1) * survival_cohort(k,i);
      cumulative_nat_mort_cohort(k) += cumulative_survival_cohort(k,i-1)*(1. - survival_cohort(k,i))*M/Z_cohort(k,i);
    }
 }
  ll.initialize();
  for(i = 1; i<=fishery_n_obs; i++)
  {
	  
	  log_exploit = 0.0;
	  ll_i = 0.0;
		  //cout << "ll_i: " << ll_i << endl;
	  if(fishdat_cohort_intrvl_index(i) > 1)
	  	ll_i += log(cumulative_survival_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)-1));
	  
	  log_exploit =  log(1-survival_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)));
	  log_exploit += log(F_fishery_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)));
	  log_exploit -= log(Z_cohort(fishdat_cohort_index(i),fishdat_cohort_intrvl_index(i)));

	  ll_i += log_exploit;
	  ll_i *= fishdat_catch(i);
      ll_i -= gammln(fishdat_catch(i)+1.);
	  ll += ll_i;
	  
  }
  for(i = 1; i<=survey_n_obs; i++)
  {
	  
	  log_exploit = 0.0;
	  ll_i = 0.0;
	  if(survey_cohort_intrvl_index(i) > 1)
	  	ll_i += log(cumulative_survival_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)-1));
	  
	  log_exploit =  log(1-survival_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)));
	  log_exploit += log(F_survey_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)));
	  log_exploit -= log(Z_cohort(survey_cohort_index(i),survey_cohort_intrvl_index(i)));
	  
	  ll_i += log_exploit;
	  ll_i *= survey_catch(i);
      ll_i -= gammln(survey_catch(i)+1.);
	  ll += ll_i;
  }

  N0 = mfexp(log_N0);
  for (k=1; k<=n_cohorts; k++)
  {
  	ll += gammln(total_fishery_catch_cohort(k) + N0(k) + total_survey_catch_cohort(k) + 1.) - gammln(N0(k)+1.);
  	ll += N0(k)*log(cumulative_survival_cohort(k,n_intrvls_cohort(k)) + cumulative_nat_mort_cohort(k));
  }	
  cout << "ll: " << ll << endl;
  nll = -1.*ll;

REPORT_SECTION
  report << "q" << endl;
  report << q << endl;
  report << "survey_sel50" << endl;
  report << survey_sel50 << endl;
  report << "survey_selslope" << endl;
  report << survey_selslope<< endl;
  report << "fishery_sel50" << endl;
  report << fishery_sel50 << endl;
  report << "fishery_selslope" << endl;
  report << fishery_selslope<< endl;
  report << "M" << endl;
  report << M << endl;
  report <<"N0"<<endl;
  report << N0 + total_fishery_catch_cohort + total_survey_catch_cohort <<endl;
  report <<"F"<<endl;
  report << F <<endl;

TOP_OF_MAIN_SECTION
  gradient_structure::set_MAX_NVAR_OFFSET(1600);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(200000);
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(1800); 
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2000000);
  arrmblsize=100000000;
